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Abstract: Considering the significance of low-energy electrons (LEEs; 0-20 eV) in radiobiology, the 
sensitization potential of gold nanoparticles (AuNPs) as high-flux LEE emitters when irradiated with 
sub-keV electrons has been suggested. In this study, a track-structure Monte Carlo simulation code 
using the dielectric theory was developed to simulate the transport of electrons below 50 keV in gold. 
In this code, modifications, particularly for elastic scattering, are implemented for a more precise 
description of the LEE emission in secondary electron emission. This code was validated using the 
secondary electron yield and backscattering coefficient. To ensure dosimetry accuracy, we further 
verified the code for energy deposition calculations using the Monte Carlo toolkit, Geant4. The 
development of this code provides a basis for future studies regarding the role of AuNPs in targeted 
radionuclide radiotherapy. 
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1 Introduction: 

The higher radiobiological significance of medium-energy (< 100 keV) and low-energy (< 20 eV) 
electrons (LEEs) compared to the commonly used high-energy (> 100 keV) electrons has been 
recently emphasized in radiotherapy. These medium-energy electrons, particularly the sub-keV Auger 
electrons (20-500 eV), are deposited across a subcellular (< 0.5 um) range in tissue!!! and exhibit a 
higher linear energy transfer (LET: 4-26 keV/uml”!) and larger relative biological effectiveness (RBE: 
1-206) than commonly used ~MeV X/y-rays (LET: 0.2-2 keV/um, RBE: 1) in radiotherapy. Thus, 
similar to heavy-charged particles, this type of electron can induce extreme toxicity in small tumors, 
clusters of cells, or micrometastases. They also have significant potential for highly targeted implant 
radiation therapy'!; however, this is only possible if they are located sufficiently close to the cell 


nucleus or membrane. More importantly, these nearly non-ionizing LEEs are generally considered 


less important in radiology owing to the ionization potential (11.7 eV!) of liquid water. Following 
the pioneering work by Sanche et al.'°!, several studies!” 8! have shown that ultralow-energy electrons 
(LEEs or quasi-free electrons (egf ): 0-20 eV) can increase the yield of DNA strand breaks. Although 
certain LEEs are unable to ionize molecules, they can still cause DNA cleavage, such as single-strand 
breaks (SSBs) and double-strand breaks (DSBs) with a high cross section, via dissociative electron 
attachment"!, Therefore, a precise description of medium-energy electrons and LEE transport in 
certain biological systems is critical. 

In targeted radionuclide therapy, radiolabeled high-Z nanoparticles, particularly gold 
nanoparticles (AuNPs), increase the absorbed dose in the near region and enhance the emission of 
LEEs through interactions with 10-200 keV photons and ~MeV electrons!!°'3!, Unfortunately, 
external beams in radiotherapy, such as photons and electrons, do not penetrate biological tissues 
beyond the mm and cm scales, respectively!!® 14+ !51, This limits the applicability of AuNPs in localized 
nanoscale radiotherapy. To address this problem, an Auger electron source attached to the AuNPs has 
been suggested!!®! to produce a larger quantity of secondary LEEs from their interactions with gold. 
These LEEs would have shorter ranges than the initial Auger electrons and further localize radiation 
energy within approximately 1-10 nml!” !81, Nevertheless, the influence of AuNPs irradiated with 
these medium- and low-energy Auger electrons has been rarely studied, and the accuracy of such 
studies cannot be guaranteed. A recent study!!*! found that a stable monolayer film of 1I on gold- 
coated mica substrates amplifies the LEE emission by six times. Owing to the abundant Auger 
electrons emitted from the decay of '7°I (23.0 Auger electrons per decay with an average energy of 0.5 
keV"°!) bombarding the gold surface, the increase in LEEs was attributed to the secondary electron 
emission (SEE) from the metal. Previous SEE experiments with a gold foil have shown that primary 
electrons of 500-750 eV generate the highest secondary electron yield (SEY; the number of true 
secondary electrons (0-50 eV) emitted per incident electron) of approximately 1.2—1.77!73!, Based 
on these experiments, the potential of LEE therapy combining Auger electron emitters with AuNPs 
has been suggested!!*!, Therefore, the primary aim of this study is to build a Monte Carlo simulation 
code for the transport of sub-keV electrons and LEEs in gold. 

In conventional Monte Carlo simulations, electron transportation is processed through the so- 


called condensed history method, including multiple scattering, single scattering, and continuous 


slowing down approximation (CSDA) strategies“. These strategies are effective and fast in 
macroscopic situations, yet not detailed in describing low-energy electrons owing to the neglect of the 
actual electron energy loss in inelastic collisions. More accurate techniques such as energy straggling 
strategy!‘ in material science or track structure strategy'*>! for microdosimetry have been developed. 
These methods require a detailed description of each electron interaction and considers the discrete 
energy loss during electron transport. There have been various well-developed generic Monte Carlo 
toolkits, such as Geant4°!, MCNP?7!, IMCT!®!, SuperMC?"!, and EGSnrc°!, which are suitable for 
the simulation of many particles. However, most generic Monte Carlo codes use the condensed history 
method rather than the track-structure method. Recently, track-structure tools or codes, such as 
Geant4-DNA?!4], NOREC®*!, KURBUC!*!, and PARTRAC7! have been implemented to simulate 
discrete electron transport in liquid water down to several eV. Sakata et al. developed two discrete 
models called Geant4 DNA Au 2016*! and Geant4 DNA Au 2018'°! for electron transport in 
gold down to 10 eV. An extension of Geant4 MicroElec was also recently developed for 
microdosimetry and SEE simulations in various materials*°!, He et al.'4!! developed a Monte Carlo 
code for low-energy incident electron transport (down to 50 eV) in Al, Ti, Pt, or Au and adopted two 
strategies, namely, the CSDA and track structure. Nevertheless, to study LEE emission from the AuNP 
surface, Monte Carlo codes considering both LEE emission from the gold surface and the dosimetry 
accuracy in gold material are essential but remain lacking. 

Accordingly, a refined Monte Carlo simulation code for the sub-keV electron transport in 
nanoscale gold material was developed. This study aims to simulate LEE emission from gold material 
when irradiated with sub-keV electrons. This would be beneficial for future research regarding the use 
of AuNPs in targeted radionuclide therapy. In this study, we demonstrate the simulation structure and 
discuss the physical mechanism and refinements in detail. A truncation method of the elastic cross 
sections is implemented by considering the Mott theory and electron-acoustic phonon scattering 
theory. In the inelastic calculations, the exchange effect, Born correction, and relativistic effects are 
simultaneously considered. The initial energy of the Fermi electrons in the conduction band electron 
excitation was also refined from the Fermi energy to the sampled energy from the joint density of 
states. The simulated SEE and calculated energy deposition were used to evaluate the reliability of the 


developed code for low-energy levels and the accuracy of the LEE emission and dosimetry 


calculations in future studies. 
2 Experimental 

The present Monte Carlo code was built for low-energy electron simulations in gold using the 
energy straggling strategy!*4!, A comprehensive process for the strategy is described by Tang et al.“”, 
This method considers the specific electron interactions and individual energy loss in each inelastic 
event, which increases the accuracy of the secondary electron yield and energy distribution calculation. 
The low energy limit of the developed code was 0.5 eV above the Fermi level (5.53 eV) of gold. 
Moreover, bremsstrahlung was not considered, and the reasonably high-energy limit of the code is 
100 keV. The basic structure of the simulation is illustrated in Fig. la. All formulae below are 
expressed in atomic units. 
2.1 Elastic Scattering 

The elastic scattering between the incident electrons and the strong nuclear Coulomb field 
represents the main cause of the angular deflection in the electron trajectory. Owing to the large mass 
of the nucleus compared to that of the electrons, the energy transfer of electrons during elastic 
scattering is negligible. In this study, the electron elastic total cross sections and differential cross 
sections were calculated based on the Mott theory'**! using the ELSEPA code developed by Salvat et 
al.“* 451, The muffin-tin radius (1.44 A) of gold was also considered in the calculation for solid gold. 
In the latest version of ELSEPA™*!, this low-energy limit was extended down to 5 eV. The modification 


of elastic scattering was introduced by the following electron-acoustic phonon scattering. 
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Fig. 1 a Schematic diagram of the simulation structure. b Phonon dispersion relations (longitudinal 
and transverse mode) for Au (circles represent experimental data)!**!, e New elastic mean free path 
(Mott-eph) interpolated from the Mott cross sections and electron-phonon cross sections 

At very low energies, electrons tend to behave like Bloch waves rather than as particles!*”!, Thus, 
they can interact with the collective movement of atoms in a solid, that is, phonons. We applied the 
theory for electron-acoustic electron scattering described by Fitting et al.'*! and modified by 


Verduin'*?! and Theulings'*"!. The electron-acoustic phonon scattering rate is given as follows: 
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where Eac is the acoustic deformation potential, kg is the Boltzmann constant, T is the 
temperature, c, is the sound velocity, pm is the mass density, and A = 5Epz is the screening factor 


adopted by Bradford and Woolf®"!, The variable Egz is the Brillouin zone energy given by Egz = 


2 
ae where the boundary wave vector is found at kgz = 2m/a if the first Brillouin zone is 
assumed to be spherical. On the other hand, a is the lattice constant, and D(E) is the density of 
states (DOS), which is typically assumed to have a parabolic dispersion relation. In this case, D(E) 


is expressed as follows: 


y 2m (E — Ecg) 
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where mp is the effective DOS mass and Ec, is the energy at the bottom of the conduction band 
and is equal to 0 for metals. 
Nez is the acoustic phonon population N(k,T) at the Brillouin zone boundary (k = kgz). The 


acoustic phonon population N(k,T) follows the Bose-Einstein distribution. Thus, 
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N(k,T) = 


where Wgz is the acoustic phonon energy W,-(k) at the Brillouin-zone boundary. The acoustic 
phonon energy w,-(k) follows a dispersion relation that was modified to fit the experimental data!**! 


as follows: 


Wac(k) = csk — ak? — BR?, (4) 


where a and p are the fitting parameters. There are three acoustic phonon modes: one longitudinal 
and two transverse. For the longitudinal mode, the chosen parameters are c, = 3240m/s, a = 
3.5x10° m? /s,and B = 3.5x 10! m3 /s. On the other hand, the chosen parameters for the transverse 
modes are c, = 1200 m/s, æ = -7.6x10° m?/s, and B = 5.3x10!8m3/s. The fitted phonon 
dispersion relations are demonstrated in Fig. 1b. 

Subsequently, the inverse mean free path (MFP) or macroscopic cross section can be derived by 


dividing the scattering rate by the electron velocity, which can be expressed as 
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In region [Egz/4,Egz|, the inverse MFP is interpolated between both values at the upper limit of 
the top formula and the lower limit of the bottom formula. According to the angular differential 


scattering rate, the angular differential inverse MFP is given by 
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The angular distribution of electron-acoustic scattering can also be obtained from the equation 
above. 

Both acoustic phonon emission and absorption exist in electron-acoustic scattering; however, the 
probability of the former is larger than that of the latter. Therefore, the weighted average energy loss 
is primarily considered. Although the energy loss is only a few meV, it cannot be neglected in low- 
energy regions. Thus, the average energy loss Eac during electron-acoustic phonon scattering can be 
approximated by the following expression: 

E fo?” ħwac(k)k?dk 
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Finally, the acoustic deformation potential of the metal was derived using the resistivity based on 


the model of Kieft and Bosch'”!. The acoustic deformation potential can then be derived as follows: 


(8) 


where pp is the resistivity. Using the parameters indicated above, the acoustic deformation potential 


of gold was found to be 4.86 eV and 1.80 eV for the longitudinal and transverse modes, respectively. 

As suggested by Fitting et al.!4’->7!, electrons should behave like Bloch waves and interact with 
phonons at low energies (< 100 eV); however, Mott cross sections are no longer valid because the 
unrealistically high scattering rate leads to a MFP that is shorter than the interatomic distance. 
Therefore, we modified the elastic scattering as follows: Mott cross sections were applied for energies 
above 100 eV, whereas electron-acoustic phonon cross sections were used for energies below 50 eV. 
The elastic cross sections between 50 eV and 100 eV were derived through interpolation between the 
two cross sections. As a result, a new elastic MFP is demonstrated in Fig. 1c. Although the elastic 
scattering models may not be sufficient for describing the structure of small AuNPs, a more precise 
atomic potential for the Mott theory will be obtained in the future using more advanced techniques, 
such as ab initio and DFT. 
2.2 Inelastic Scattering 

This study mainly considers electron-electron inelastic scattering. Inelastic cross sections were 
computed using the dielectric theory'>*!. The core of this theory is the complex dielectric function 


e(q, w), where q is the momentum loss and w is the energy loss. The imaginary part of the inverse 
dielectric function is called the energy loss function (ELF) Im Econ al which provides detailed energy 


loss information for the material electrons. From the optical data of experiments!*“! or ab initio and 
DFT calculations", the dielectric function and ELF can only be obtained at the optical limit (q > 0), 
which is called the optical dielectric function and optical ELF. 

To extend the ELF from the optical limit to q > 0, we used the extended Drude model!**! and 
Mermin model!°’!, In these two models, the sum of the Drude-type ELFs in (4) is used to fit the optical 
ELF because both the extended Drude model and Mermin model coincide at the optical limit!®!, that 


is, the Drude-type ELFs: 
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where A; = a;w2, is the oscillator intensity, y; is the damping constant, and Wp, is the plasma 


p.i 


frequency. 

The fitting model was built based on several previous studies!°* 4% 5° ©], In the development of 
Geant4 MircoElec, Valentin'°! and Gibaru'*®! et al. used a Heaviside step function to cut each 
oscillator in the fitting. Pablo et al.'°"! used an exponential smooth switching function instead of a 
simple step function for the outer-shell electrons. Based on these prior methodologies, the fitting 


model is expressed as follows: 
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For the outer-shell electrons, F(w) represents a smooth exponential switching function for the 
electronic excitations!*"): 
1 
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where A; and Wp; are the additional fitting parameters. For the inner shells, particularly the K and 


L shells, which have sharp edges in the optical ELF, F(w) is chosen as a simple step function 


O(w = wth i). For plasmons and conduction band electrons, F (w) can be simply set to 1. 
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Fig. 2 Comparison of the fitted ELF (solid) and experimental ELF (circles): a contribution of each 
electronic shell obtained from the combined oscillators at the energy loss range of 0-100 eV; b 
contribution to all electronic shells of the combined oscillators 

The experimental optical ELF, which covers the electron energies from 0.005 eV to 1 MeV, is 
derived from the databases of Palik!®?!, Windt'®?!, and the evaluated photo data library 1997 version 


(EPDL97'™!), in which the high-energy optical ELF is calculated using the photo-ionization cross 


section!®! as follows: 
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where nm is the atomic density, c is the speed of light, and Gpnot is the photoionization cross 
section. 

The fitting results were validated by using two sum rules (f-sum and ps-sum) and the mean 
ionization potential. A partial sum rule for the fitted oscillators was also used to distinguish the 
contribution of the oscillators to the electronic shells. Initial values of the fitting (energy loss from 0 
to 100 eV) were partially obtained by referring to the study of Pablo et al.'*!!. We managed to use 28 
oscillators to describe the optical ELF. Fig. 2 demonstrates that the fitted ELF data sufficiently agrees 
with the experimental ELF. 

During fitting, the plasmon peak was chosen as 25.42 eV, which is near 24.8 eV in Egerton’s 
book!**!, We verified the sum rule of the fitted ELF, that is, Zef = 78.05 (ideal value: 79.0), while 
the ps-sum rule resulted in 1.0908 (ideal value: 1.0). The mean ionization potential based on the fitted 
ELF was 736.45 eV. These values obtained from the fitted ELF coincide with those from the 
experimental ELF, with only a small deviation in the high-energy region (L and K shells). This is 
likely due to the restriction of the Drude functions. If a more accurate experimental ELF is available, 
then better fitting results can be obtained. We also checked the partial sum rule to determine the shell 
to which the fitted oscillators belong, and the results are also displayed in Fig. 2. Therefore, these 
categorized oscillators enabled the calculation of the inelastic mean free path (IMFP) of the different 
shells. 

The calculation of IMFP is achieved by extrapolating the momentum loss q from the optical limit 
to q > 0. For the extended Drude model, the extrapolation of g was achieved by using the dispersion 
relation. The extrapolation was performed using the following equation: 
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In this method, w,j;(q) describes the plasmon dispersion and y;(q) is the g-dependent damping 

constant, which represents the broadening effect. The variable y; is usually considered a constant; 


6Ef 


2 2 
thus, one can simply choose (1) @yi(q) = @pi + T or (2) Wpi(q) = opi tq + L [56,67]. If 


Yi is considered a function of q, another dispersion relation is given by! ®! (3) wyi(q) = 


Ws i + r q2 + Z and y;(q) = lv? + £, This was already presented in the study of Ritchie et 
al.°*! and proved!*! to be more comprehensive than the extrapolation considering only the w, ;(q). 
Therefore, we chose the third extrapolation method for the conduction band and outer shells and the 
first method for the inner shells as the first inelastic model, which is labeled as Ritchie in the following 
text. 

Regarding the Mermin model, the extrapolation procedure was already achieved using the Mermin 
dielectric function!””! as follows: 

(14+) [eq +iy)-1 
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where ¢, is the Lindhard dielectric function!’ and ¢,(q,0) is the Lindhard dielectric function at 


Eu(q,w) =1+ (14) 


the static limit. This is the second inelastic model, which is also labeled as Mermin. 

A significant advantage of the extended Drude model and Mermin model is that fitted oscillators 
can approximate the contribution to particular electronic shells or the plasmon in the calculation. 
Therefore, combined with the divided shells in the ELF, we can calculate the IMFP of the different 
electronic shells. 


The general calculation of IMFP considering the relativistic high energy correction is performed by 
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where E is the kinetic energy of the incident electrons, c is the speed of light, E is the 
2c2 


relativistic correction, W, is the maximum energy loss, W_ is the minimum energy loss, and q4 


and q_ are the maximum and minimum momentum losses, respectively. The variables q, and q_ 


q4 = e(2+5)4 e—a)(2+#59), (16) 


2 2\2 
where fex(q,E) =1—- c + (5) is the Ochkur exchange correction”. For a collective excitation, 


are given by 


fex = 1. As indicated above, F(w) is an exponential switching function or a simple step function. 
The inelastic cross sections of each shell were individually calculated. For metals such as Au, W_ 


equals zero. As proposed by Pablo et al.!6'], the maximum energy loss W, depends on the different 


electronic shells. Plasmon and 5d electrons are considered collective excitations; thus, W} = E — Ep. 
For the outer-shell and inner-shell electrons, W} equals min[E — Ep, (E + B;)/2], where B; is the 
binding energy of the i-th electronic shell. The binding energies of the outer- and inner-shell electrons 
were derived from the FFAST database!”*!, Meanwhile, a classical Coulomb-field Born correction that 
accounts for the potential energy gained by the incident electrons is added for the calculation of the 
outer-shell and inner-shell electrons. This is achieved by replacing the kinetic energy E with E’ = 
E + 2B; in the calculation; however, the energy loss should remain unchanged. 

Based on the aforementioned discussion, we compared calculated IMFP with previously reported 


experimental!’*’°! and theoretical data!’”**!, which are shown in Fig. 3. 
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Fig. 3 a Calculated inelastic mean free path and comparison with other published data!’*”°!, b 
Calculated inner shell cross sections (K, L, M shells); solid, dashed, and dotted lines indicate the 
Ritchie model, Mermin model, and Gryzinski theory, respectively!*” 83] 

Fig. 3a presents the calculated IMFP, which sufficiently agrees with the experimental and previous 
theoretical data. A relativistic effect is clearly observed at high energies. At low energies below 100 
eV, the modified Ritchie model can derive a good IMFP similar to that of the Mermin model due to 
its improved estimation of the plasmon lifetime by considering y as a function of q. We also 
calculated the inner-shell ionization cross sections and compared them with the relativistic Gryzinsk1’s 


al87. 88] for validation. As shown in Fig. 3b, the calculation using the Ritchie model can yield 


formul 
better inner-shell ionization data compared to the Mermin model, which indicates that the Ritchie 
model is recommended for high-energy electron simulations. This may be due to the fact that for the 


inner-shells, the plasmonic consideration in the Mermin model is no longer valid, whereas in the 


Ritchie model, the simple extrapolation of q without considering the plasmon is more applicable. 
Although the inner shell cross sections using both the Ritchie and Mermin models remain lower than 
those calculated by the Gryzinski theory, which is considered more accurate, this may be due to the 
fact that the fitted ELF is insufficiently accurate for the inner shells. Therefore, this lowers the 
accuracy of the calculated inner-shell cross sections. 

Following an inelastic event, in normal considerations, a secondary electron will be generated with 
an energy of Es = W + Ep for Fermi electron excitation (conduction band and plasmon excitation as 
considered in this study) and Eş = W — B; for the inner- and outer-shell ionization, where W is the 
energy loss. In this case, we applied the improvement made by Ding et al.'**!, which assumes that the 
initial energy E’ of the Fermi electron is proportional to the joint density of states (JDOS) of the free- 
electron gas, that is, P(E’, W) « JE'(E'+W) . In the expression, E’ must satisfy E’ < Ep and 
E' + W > Ep. Therefore, the modified energy of the Fermi electron after inelastic scattering is Es = 
WHE". 

2.3 Surface Potential 

For an electron to escape from the surface of metals, it must have a kinetic energy larger than the 
Fermi level Er and the work function Wp, that is, the potential barrier between the metal surface and 
vacuum"!, Therefore, the cut-off energy in the simulation is set to Ep + Wp, and the electrons are 
locally terminated below this limit, which can also save memory. 

To escape from the surface, once an electron reaches the gold surface, its kinetic energy must satisfy 


the following condition: 
Ecos*@ >U, (17) 


where E is the energy that considers the Fermi level as the zero point, 0 is the angle between the 
electron motion direction and is normal to the surface before its emission, and U = Wp is the 
potential barrier above the Fermi level. Subsequently, according to the quantum theory, the emission 


89-91] which is given by the following: 


44/1 — U/(Ecos?6) 


[1 +y1-— U/(Ecos*6)| ) 


will depend on the transmission coefficient T! 
(18) 


Arandom number R is generated for comparison with T to determine the emission as follows: 


eee if R<T (19) 


Reflection, else ` 


The emitted electrons have a kinetic energy of E — U, and electrons entering the solid will also 
gain potential energy Wp above the Fermi level. In the SEE simulation, the work function between 
the gold surface and vacuum was set to 5.1 eV. Considering practical cases of water-gold interfaces 
for future development, we will use the difference between the work functions of these two materials 
and assume that the calculation of the transmission coefficient is similar to that considered in the 
surface-vacuum interface. 

3 Results and Discussion: 
3.1 Secondary electron emission 
3.1.1 Secondary electron yield 

To validate the code, we simulated SEE experiments and calculated the secondary electron yield 
(SEY) and backscattering coefficient (BSC) of gold. This simulation was simply set up as multiple 
groups of mono-energetic electron beams vertically bombarding an infinite gold plane, and the 
information of the electrons emitted from the surface was collected. The SEY is defined as the ratio 
of the number of emitted secondary electrons (< 50 eV) to the number of incident primary electrons. 
The BSC is defined as the ratio of the number of backscattering electrons (> 50 eV) to the number of 
incident primary electrons. To obtain the SEY and BSC, 1x106 particles were simulated. The 
uncertainty of this code in the SEE simulation is less than 1%. As presented in Fig. 4, we demonstrate 
the calculated SEY and BSC and compare them with the experimental data collected by Joy et al.!°7!, 
Regarding the experimental data, monoenergetic electron beams were used in most of the experiments 


to bombard the material under vacuum conditions, which was similar to the simulation setup. 
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Fig. 4 a SEY and BSC derived from the SEE simulations and compared to the experimental data!””!, 
M indicates the Mermin model and R indicates the Ritchie model. b Comparison of the SEY values 
calculated by two ways (solid lines: Mott-eph, dashed lines: Mott) 

As shown in Fig. 4a, the simulated SEY values mostly coincide with the experimental values!””!, 
which is essential for the calculation of the emission of LEEs from AuNPs. The Ritchie model derived 
a larger SEY compared to the Mermin model, whereas the BSC values were similar in both models. 
In addition, if the modification for low-energy elastic scattering was not added, the SEY values 
calculated in both ways would mainly differ below 1 keV, as shown in Fig. 4b. The figure shows that 
the calculated SEY using only the Mott cross sections is larger than the experimental values for the 
primary electrons with a lower energy (< 400 eV). The modification of the elastic cross sections 
improved the description of the SEE. The high SEY also suggests that the '”°I emitted electrons can 
enhance the LEE emission from gold nanoparticles. 

3.1.2 Secondary electron spectra 

The secondary electron spectra are critical for the analysis of LEEs emitted from AuNPs because 
the cross sections of the LEEs damaging the DNA, such as DSBs and SSBs, are highly energy- 
dependent!?!, The 1x10” electrons were simulated to obtain stable secondary electron spectra, which 
are presented in Fig. 5. The figure demonstrates that most of the emitted secondary electrons were 


below 30 eV; a zero-loss elastic peak is also observed herein. 
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Fig. 5 Emitted electron spectra from SEE simulations: a 200 eV primary electrons, b 500 eV 
primary electrons, ¢ 1 keV primary electrons, and d 2 keV primary electrons 
Moreover, the peak of LEEs was found to be located at 2.25—3.25 eV, which is similar to the 2.2 eV 
peak value reported in the ion-bombardment SEE experiment'”*!, This large amount of LEE emission 
also suggests the feasibility of using sub-keV electrons combined with AuNPs as a high-flux LEE 


source. 
3.2 Dosimetry test 


Dosimetry testing is critical for the evaluation of further microdosimetry. The stopping power, S(E), 
for electrons of energy E can be calculated using a method similar to the IMFP. The expression for 


S(E) is shown below: 


EY : 
S(E) = oe odo f fax (q, mirom =]. (20) 
2c2 T 


As shown in Fig. 6, we compared the calculated stopping power with those of other theoretical 
calculations! } 85 94971 and experiments!’®!. The effect of the relativistic corrections introduced in the 
calculation is demonstrated by the stopping powers at high electron energies, which is in good 
agreement with the ESTAR'®*! and the evaluated electron data library (EEDL) databases!*"!, With the 
modifications, the Ritchie model can produce an accurate stopping power comparable to that of the 


Mermin model at low energies. 
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Fig. 6 Calculated stopping power and comparison with other published data!*! 85- 94-971 


The energy deposition was mainly tested using the semi-infinite gold model. Calculation of the 
energy deposition depth profiles can provide information regarding the energy deposited in gold and 
its uniformity. This verification also guarantees the correctness of dosimetry for future studies. 
Electron energy deposition as a function of the depth of semi-infinite gold was calculated for 
comparison using the developed code and Geant4 standard electromagnetic physics (option 4). In the 
calculation, 200 eV, 500 eV, 1 keV, 2 keV, and 10 keV electrons were simulated; 5x10° electrons were 


adopted in each test. The calculated results are shown in Fig. 7. 
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Fig. 7 a Calculated total energy deposition. b Calculated energy deposition as a function of depth (1 
— 1000 A with 1 A as the thickness of each layer) using the developed code (solid and dashed lines 
indicate the Mermin and Ritchie models, respectively) and Geant4 (dotted line) 

The outcomes of both inelastic models in the code present a similar tendency. Fig. 7 (a) 
demonstrates that the ratio of the total energy deposition results of this code to that of Geant4 was 
approximately 89%—93%. This slight decrease was due to the larger secondary electron emission from 
the surface. In contrast, the difference in the energy deposition depth profiles is shown in Fig. 7 (b). 
At low energies, our inelastic cross sections differ from those of the Livermore model (EEDL!*!)). 
Gibaru et al.“ suggested that the Geant4 standard physics underestimates the range of electrons, 
which is also observed upon comparing our results. At high energies, the energy deposition depth 
profiles of both this code and Geant4 tended to coincide. Note that the more detailed energy-straggling 
strategy in this code is significantly slower than the multiple scattering strategy in Geant4 standard 
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A function for calculating the energy deposition in the 3D voxels was also added to enable future 
microdosimetry studies. The method for calculating the energy deposition in the 3D voxels is 
demonstrated, and the calculation is displayed in Fig. 8. The length of the 3D voxels was set to 1 A, 


and the maximum length of each axis was 1000 A. The 1x10’ electrons were simulated here. 
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Fig. 8 a Diagram of the calculation of the energy deposition in 3d voxels; projection of the 
calculated energy deposition (Ritchie) by the 3D voxels. b 1 keV X-Z plane, ¢ 2 keV X-Z plane, d 
10 keV X-Z plane, e 1 keV X-Y plane, f 2 keV X-Y plane, and g 10 keV X-Y plane 

In Fig. 8, the spatial distribution of the energy deposition demonstrates good symmetry, which 
verifies the particle transport in the Monte Carlo procedure. The intensity of the energy deposition in 
the X-Z plane is more comprehensive than that of the depth profiles. At low energies, the projection 


of energy deposition in the X-Z plane is an ellipse around the entry point of the electron beams. As 


the incident electron energy increases, the energy deposition distribution tends to become more 


cylindrical. The distribution of the energy deposition also demonstrates the nanoscale precision of 
this code. This function allowed the further visualization and analysis of the electron energy 


deposition inside the AuNPs when irradiated with electrons. 
4 Conclusion 


A Monte Carlo code was built to simulate the nanoscale electron transport in gold material down 
to 0.5 eV above the Fermi level. Mott elastic scattering cross sections were calculated using the 
ELSEPA code!4!, Electron-acoustic phonon cross sections were calculated from the fitted dispersion 
relation. The elastic cross sections were modified using a truncation method that applies Mott cross 
sections at energies higher than 100 eV, electron-phonon cross sections at 0—50 eV, and interpolation 
in between. The optical ELF data were fitted into 28 Drude oscillators, which were categorized into 
different electronic shells to calculate the cross sections using the dielectric theory. Two inelastic 
models, the Ritchie and Mermin models, were used to produce inelastic cross sections. The exchange 
effect, Born correction, and relativistic effects were carefully considered. The initial energy of the 
Fermi electrons in the conduction band electron excitation was refined from the Fermi energy to the 
sampled energy from the JDOS. 

We simulated and compared the SEE results with the experimental results and found a reasonable 
agreement. The quantum model of the surface potential plays an important role in the calculation of 
the SEE from the surface. The simulated secondary electron spectra demonstrate a high yield of LEEs 
from the gold surface when irradiated with keV electrons. Thus, building such a simulation code is 
significant for studying LEEs in micro- and nanodosimetry. The inelastic stopping power and energy 
deposition depth profiles were validated. The 3D voxel function was implemented in the developed 
code to allow further usage. 

This code is planned for the evaluation of the LEE emission in the '*°I/AuNPs system. It is also 
being developed for its significant potential in future microdosimetry studies. Further improvements 
in this code include considering the structure of small AuNPs, implementing the surface plasmon 
mode for certain AuNPs, adding more materials, and using an advanced geometry intersection 
algorithm. 
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